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ABSTRACT 

In this work we introduce a new family of ten-step linear multistep methods for the 
integration of orbital problems. The new methods are constructed by adopting a new 
methodology which improves the phase lag characteristics by vanishing both the phase 
lag function and its first derivatives at a specific frequency. The efficiency of the new 
family of methods is proved via error analysis and numerical applications. 
PACS: 02.60, 02.70.Bf, 95.10.Ce, 95.10.Eg, 95.75.Pq 

Subject headings: Numerical Solution, N-body problem, Multistep Methods, Phase Lag, 
Phase Fitting 



1. Introduction 

The numerical integration of systems of ordinary differential equations with oscillatory solu- 
tions has been the subject of research during the past decades. This type of ODEs is often met in 
real problems, like the N-body problem. For highly oscillatory problems standard non-specialized 
methods can require a huge number of steps to track the oscillations. One way to obtain a more 
efficient integration process is to construct numerical methods with an incre ased algebraic o rder, 



although the implementation of high algebraic order meets several difficulties iQuinlanl (| 19991 ) . 



On the other hand, there are some special techniques for optimizing numerical methods. 
Trigonometrical fitting and phase-fitting are some of them, producing methods with variable coeffi- 
cients, which depend on v = coh, where oj is the dominant frequency of the problem and h is the step 
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length of integration. More precisely, the coefficients of a general linear method are found from the 
requirement that it integrates exactly powers up to degree p + 1. For problems having oscillatory 
solutions, more efficient methods are obtained when they are exact for every linear combination of 
functions from the reference set 



{l,x,...,x K ,e ± ^,...,x p e ± ^ x } 



(1) 



This technique is known as ex ponential (or trigonometric if /U = iuj) fitting and has a long history 
Gautschil (|196ll ). iLvchd (|1972l ). The set (TjQ) is characterized by two integer parameters, K and P . 
The set in which there is no classical component is identified by K = — 1 while the set in which there 
is no exponential fitting component (the classical case) is identified by P = —1. Parameter P will be 
called the level of tuning. An important property of exponential fitted algorithms is that they tend 
to the classical ones when the involved frequencies tend to zero, a fact which allows to say that ex- 
ponential fitting represents a natural extension of the classical polynomial fitting. The exami nation 
of the convergence of exponentially fitted multistep methods is included in Lyche's theory iLvche 
(|1972l ). There is a large number of significant methods pres ented with hi g h practical importanc e 



that have been presente d in the bibliography (see for example ISimo 2 d2nnd),[chawla fc Baol dl98fih 




Mazzia et al 



Raptis fe Allison ( 1978 ). Anastassi Sz Simos ( 2004 ). Anastassi Simosl ( 2005ah . lAnastassi Simos 



Anastassi fc Simosl (l2005bh.lLambert fc Watson] dl976ll.lCash & Mazzial (12 006) 



The general theory is presented in detail in llxaru Berghd ((2004). 



(|2006l 1. lBerghe Daeld (12006T) . iPsihoviosI (l2006l ). ISimosl (^OolnJSimod (|2007T ). 



Iavernaro et al 



Considering the accuracy of a method when solving oscillatory problems, it is more appropriate 
to work with the phase-lag, rather than the principal loc al truncation error. We mention the 
pioneering paper of Brusa and Nigro iBrusa Nigrd (ll980h , in which the phase-lag property was 
introduced. This is actually another type of a truncation error, i.e. the angle between the analytical 
solution and the numerical solution. On the other hand, exponential fitting is accurate only when 
a good estimate of the dominant frequency of the solution is known in advance. This means that 
in practice, if a small change in the dominant frequency is introduced, the efficiency of the method 
can be dramatically altered. It is well known, that for equations similar to the harmonic oscillator, 
the most efficient exponential fitted methods are those with the highest tuning level. 

In this paper we pre sent a new family of m e thods based on the 10-step linear multistep method 
of Quinlan and Tremain iQuinlan Tremaind ([1990 ). The new methods are constructed by van- 
ishing the phase-lag function and its first derivatives at a predefined frequency. Error analysis and 
numerical experiments show that the new methods exhibit improved characteristics concerning the 
long term behavior of the solution in the 5-outer problem. The paper is organized as follows: In 
section 2, the general theory of the new methodology is presented. In section 3, the new methods 
are described in detail. In section 5 the stability properties of the new methods are investigated. 
Section 5 presents the results from the numerical experiments and finally, conclusions are drawn in 
section 6. 
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2. Phase-lag analysis of symmetric multistep methods 



Consider the differential equations 

= f(t, y), y(to) = y , y'(t ) = y' Q (2) 

and the linear multistep methods 

J J 

^2 ajVn+j = h 2 ^ b jfn+j (3) 
3=0 3=0 

where y n+j = y(t + (n + j)h), f n+j = f(t + (n + j)h, y(t + (n + j)h)) and h is the step size of 
the method. We associate the following functional to method ([3]): 

J J 
L(h, a, b, y{t)) = ajy(t + j-h)-h 2 ^2 &*!/"(* + 3 ' h) (4) 

3=0 J=0 

where a, b are the vectors of coefficients aj and bj respectively, and y(t) is an arbitrary function. 
The algebraic order of the method ([3]) is p, if 

L(h,a,b,y(t)) =C p+2 h p+2 y^{t) + 0{h^) (5) 

The coefficients C q are given 

c o = Ej=o a 3 

C g = 7 >: ; ' u/'-" ; T^>:/ u/' % (6) 
The principal local truncation error (PLTE) is the leading term of © 

plte = C p+2 h p+2 y^ +2 \t) (7) 

The following assumptions will be considered in the rest of the paper: 

1. a j = 1 since we can always divide the coefficients of ([3]) with aj. 

2. \oq \ + \bo\ ^ since otherwise we can assume that J = J — 1. 

3. E/=o 1^3 1 ® smce otherwise the solution of © would be independent of ([2]). 

4. The method (|3|) is at least of order one. 

5. The method (J3j) is zero stable, which means that the roots of the polynomial 

J 

p(z)=J2^z j (8) 

3=0 

all lie in the unit disc, and those that lie on the unit circle have multiplicity one. 
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6. The method (J3j) is symmetric, which means that 

a j = a J-j> b j = b J-j> 3 = °( 1 ) J ( 9 ) 

It is easil y proved then that both t he order of the method and the step number J are even 
numbers lLambert fc Watsonl (|1976l ). 

Consider now the test problem 

y"(t) = -u> 2 y{t) (10) 

where a; is a constant. The numerical solution of (|1U|) by applying method ^ is described by the 
difference equation 

J/2 

M s2 )(Vn+j + Vn-j) + A (s 2 )y n = (11) 

i=i 

with 

Aj{s 2 ) = a 1 _ 1 +s 2 -bj^ (12) 

2 J 2 J 

and s = ojh. The characteristic equation is then given 

J/2 

^Ajis 2 )^ + z-i)+A (s 2 ) = (13) 
j=i 

and the interval of periodicity (0, Sq) is then defined such that for s £ (0, so) the roots of (fl~3j) are 
of the form 

Zl = e iX{s \ z 2 = e~ iX{s \ \ Zj \ < 1, 3 < j < J (14) 
where A(s) is a real function of s. The phase-lag PL of the method ([3]) is then defined 

PL = s- A(s) (15) 

and is of order q if 

PL = c • s 9+2 + 0(si +4 ) (16) 

In general, the coefficients of the method ([3]) depend on some parameter v, thus the coefficients 
Aj are functions o f both s 2 and v. The following theorem has been proved by Simos and Williams 



Simos &: Williams! (|1999l ): For the symmetric method (|10p the phase-lag is given by 

DTt , 2Z-L 2 1 Ms 2 ,v)-cos(j-s) + A (s 2 ,v) 

PL(s,v) = J - m (17) 

2E/2i 2 ^(« 2 ^) 

We are now in position to describe the new methodology. In order to efficiently integrate oscillatory 
problems, a good technique is to calculate the coefficients of the numerical method by forcing the 
phase lag to be zero at a specific frequency. But, since the appropriate frequency is problem 
dependent and in general is not always known, we may assume that we have an error in the 
frequency estimation. It would be of great importance to force the phase-lag to be less sensitive 
to this error. Thus, beyond the vanishing of the phase-lag, we also force its first derivatives to be 
zero. 
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3. Construction of the new methods 
3.1. Base Method 

The family of new methods is b ased on the 10-step linear multistep method of Quinlan and 
Tremain I Quinlan Tremaind (|1990l ) which is of the form ([3]) with coefficients 



Oo = 1 CL\ = — 1 02 = 1 «3 = — 1 O4 = 1 05 = —2 

h = n h — 399187 t _ 17327 h _ 597859 l 704183 h _ 465133 

u ° u U1 241920 " 2 8640 ° 3 60480 4 60480 " 5 24192 

The principal term of the local truncation error (PLTE) of the method is given 



(18) 



p l te= ^^y (12)fcl2 (19) 

1 912384 y v ; 



3.2. Method PF-DO: Phase fitted 

The first method of the family (PF-DO) is constructed by forcing the phase-lag function to be 
zero at the frequency V = u * h. Coefficients a are left the same, while coefficients b become: 

U n U _J l,num 1 I 2, num. 

°0 — u > u l — 8064 D ' " 2 ~ 4032 D Q ' 

1 ^3,num L 1 V™ m I* 1 ^5,num 



L J 3,num 1 . 



where 



2016 D ' " 4 4032 D ' " s 4032 D 



D = w 2 ((cos(w)) 4 + 6 (cos(w)) 2 + 1-4 cos(u ) - 4 (cos(w)) 3 ) 



and 



6° num = -16128 (cos(f )) 3 + 45139 (cos(v))V + 6048 (cos(t> )) 2 - 73215 (cos(u))V 
+3024 cos(v) + 47553 cos(v)v 2 - 11917 u 2 + 16128 (cos(v)) 5 - 8064 (cos(v)) 4 - 1008 
h \num = -32256 (cos(u)) 4 + 45139 (cos(u)) V - 64512 (cos(w)) 3 + 24192 (cos(t;)) 2 
-22026 (coe(u))V + 9656 cos(v)v 2 + 12096 cos(v) - 2529 v 2 - 4032 + 64512 (cos(v)) 5 
b\ num = 56448 (cos(v)) 4 - 73215 (coe(u))V + 112896 (cos(u )) 3 - 23113 (cos(u )fv 2 
-42336 (cos(v)) 2 + 73215 (cos(u))V - 40011 cos{v)v 2 - 21168 cos(v) + 10204 v 2 
+7056 - 112896 (cos(v)f 

b\num = -225792 (cos(u)) 4 + 325629 {cos{v)) A v 2 - 451584 (cos(u )) 3 - 38624 (cos(v)) 3 u 2 
-96246 (coe(w)) V + 169344 (cos(u)) 2 + 28968 cos{v)v 2 + 84672 cos{v) + 451584 (cos(v)) 5 
-8047 v 2 - 28224 

&° num = -388196 (cos(u))V + 282240 (cos(u)) 4 + 564480 (cos(v)) 3 - 27081 {cos{v)fv 2 
+233349 (cos(i;))V - 211680 (cos(tj)) 2 - 111571 cos{v)v 2 - 105840 cos(u) + 28899w 2 
-564480 (cos(w)) 5 + 35280 



Since for small values of v, the above formulae are subject to heavy cancelations, we give the 
Taylor expansions of b coefficients: 
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_ 399187 52559 2 100673687 4 1084493 ^ 96453547 , 

1 ~ 241920 ~ 912384 V + 29059430400 V ~ 27897053184 V + 213412456857600 V 
_ 17327 52559 2 100673687 4 1084493 6 96453547 8 

2 ~ ~ 8640 + 114048 V ~ 3632428800 V + 3487131648 V ~ 26676557107200 V + 
_ 597859 367913 ^ 100673687 ^ 1084493 ^ 96453547 8 

3 ~ 60480 ~ 228096 V + 1037836800 V ~ 996323328 V + 7621873459200 " 

_ 704183 367913 2 100673687 4 1084493 6 96453547 8 

4 ~ ~ 60480 + 114048 V ~ 518918400 V + 498161664 V ~ 3810936729600 V + " 
_ 465133 1839565 2 100673687 4 5422465 6 96453547 8 

5 ~ 24192 ~ 456192 V + 415134720 V ~ 1992646656 V + 3048749383680 V 



The PLTE of the method is given 

pUe = ( ^ (10) ^2559 (12) \ 12 

F V 912384 912384 y J y J 



3.3. Method PF-D1: Phase fitted + 1st Derivative of Phase-lag is zero 

The second method of the family (PF-D1) is constructed by forcing the phase-lag function and 
its 1st derivative to be zero at the frequency V = oj * h. The coefficients a are left the same, while 
coefficients b become: 



bo - bi - 192-^— b 2 - 

1,1 lyl ijl 

L ',_ 3, num. L I 4,num ^ l_ 5, num. 

03 ~~ 48 Di " 4 ~~ 48 Di " 5 ~~ 96 Di 

where 



£>i = t, 3 ((cos (v) f - 3 (cos (w)) 4 + 2 (cos (v)) 3 + 2 (cos (v)) 2 - 3 cos (v) + 



and 



b\ num = —48 sin(-y) — 768 (cos(u)) 3 sin(v) + 144 cos(v) sin(v) + 768 sin(u)(cos(i;)) 5 
-384 (cos(w)) 4 sin(u) + 288 (cos(v)) 2 sin(u) + 432 cos(v)v - 941 cos(v)t; 3 + 281 v 3 
+1344 (cos(v)) 5 t; + 1344 (cos(w)) 4 t; - 1776 (cos(w)) 3 w - 768 (cos(v)) 6 v - 576 (cos(v)) 2 v 
+259 (cos(w)) V + 1481 (cos(v)) 3 v 3 
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b 2,num = 72 sin ( v ) - 384 sin(i;)(cos(t;)) 6 - 960 sin(t>)(cos(t;)) 5 + 960 (cos(u)) 4 sin(u) 

+1008 (cos(v)f sin(u) - 504 (cos(t;)) 2 sin(v) - 192 cos(i>) sin(u) - 1305 (cos(u )) 3 w 3 

+732 (cos(w)) 2 w + 427 (cos(v)) V + 225 cos(v)t; 3 - 612 cos(v)v - 12 t> - 26 v 3 

+192 (cos(w)) 7 ?; + 576 (cos(w))S - 2304 (cos(v)) 5 v - 1296 (cos(v)) 4 v - 1481 v 3 (cos(v)) 4 

+2724(cos(w)) 3 t; 

& 3,n«m = - 192 sin ( v ) + 2304 sin(t>)(cos(t;)) 6 + 1920 sin(v)(cos(V)) 5 

-3840 (cos(v)) 4 sin(u) - 2208 {cos(v )) 3 sin(v) + 1584 (cos(v)) 2 sin(u) + 432 cos(v) sm(v) 

-7464 (cos{v)) 3 v + 3387 (cos(v))V - 1512 (cos(v)) 2 v - 895 (cos(v)) V 

-1088 cos{v)v 3 + 1512 cos(v)t; + 72 v - 1152 (cos(v)) 7 v + 384 (cos(u )) 6 t> + 7104 (cos{v)) 5 v 

+1481 w 3 (cos(t>)) 5 + 4443 v 3 (cos(v)) 4 + 1056 (cos(v)) 4 t> + 232 v 3 

b\ num = -192 sin(u) + 2304 sin(w)(cos(t;)) 6 + 1920 sin(t;)(cos(?;)) 5 

-3840 (cos(v)) 4 sin(u) - 2208 (cos(v)) 3 sin(t>) + 1584 (cos(v)) 2 sin(u) + 432 cos( v) sin(u) 

-7464 (cos(t;)) 3 t; + 3387 (cos(u))V - 1512 (cos(t;)) 2 t; - 895 (cos(u)) V 

-1088 cos{v)v 3 + 1512 cos(t>)t> + 72 v - 1152 (cos(v)) 7 v + 384 (cos(u )) 6 t> + 7104 (cos{v)) 5 v 

+1481 w 3 (cos(t>)) 5 + 4443 v 3 (cos (v)) 4 + 1056 (cos(t;)) 4 t; + 232 1; 3 

6| num = -720 sin(u) + 15360 sin(t;)(cos(w)) 6 + 3840 sin(u)(cos(t;)) 5 

-21120 (cos(v)) 4 sin(v) - 5760 (cos(u )) 3 sin(u) + 7200 (cos(w)) 2 sin(v) + 1200 cos(t;) sm(v) 

-29040 (cos(v)) 3 t; + 10059 (cos(v)) V - 3360 (cos(u )) 2 v - 4231 (cos(i;)) V 

+5040 cos{v)v - 3411 cos(v)w 3 + 480 v - 7680 (cos(v)) 7 t; + 11520 (cos(v)fv 

+31680 (cos(v)) 5 t; + 12252 v 3 {cos(v) f - 8640 {cos(v)) 4 v + 22484 v 3 (cos(f)) 4 + 647 v 3 

Since for small values of v, the above formulae are subject to heavy cancelations, we give the 
Taylor expansions of the coefficients b: 



399187 52559 9 

bi = f + 

241920 456192 
17327 52559 , 

b 2 = 1 v - 

8640 57024 



975124291 4 2896813 fi 1818828019 s 

v 4 v 6 H v 8 

174356582400 49816166400 1067062284288000 
M.^i 2 4461254807 , 2517959 6 98779707713 8 

2 ~ ~ 8640 + 57024 V ~ 43589145600 V + 340540200 V ~ 266765571072000 V + "' 
597859 367913 , 3127415341 4 3766196569 „ 44891085091 8 

b 3 = v 2 H v 4 v 6 H v s - ... 

60480 114048 6227020800 87178291200 20520428544000 
704183 367913 9 7330976207 d 584032469 fi 1452594367391 „ 

6 4 = 1 v 2 v 4 H v 6 v s 

60480 57024 6227020800 5448643200 266765571072000 
_ 465133 _ 1839565 2 3844845691 4 _ 4974280813 6 773868209533 8 

5 ~ 24192 ~ 228096 V + 2490808320 V ~ 34871316480 V + 106706228428800 V 
The PLTE of the method is given 

p l te = ( w 2 3,(10) + w * (8) + (12) A fc i2 (21) 

^ V 456192 912384 y 912384 y / V ; 
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3.4. Method PF-D2: Phase fitted + 1st, 2nd Derivatives of Phase-lag are zero 

The third method of the family (PF-D2) is constructed by forcing the phase-lag function and 
its 1st and 2nd derivatives to be zero at the frequency V = lo * h. The coefficients a are left the 
same, while for the coefficients b we have: 



bo 
b 3 



bi 

1 b% 

1 6,num 



bi 



b 2 

1 l.num 7, 

32 Dl "2 
h 2 

1 4,num i 

16 D 2 " 5 



h 2 

1 u 2,num 
16 D 2 

h 2 

1 5,num 

16 D' 



where 



D 2 
Do 



- „4 



- „4 



((cos (v )) 5 - 3 (cos (w)) 4 + 2 (cos (v)f + 2 (cos (v)) 2 - 3 cos (v) + l) 
(sin (v)) 4 ((cos (v)) 2 - 2 cos (v) + l) 



and 



b\ num = -48 sin(-y) - 768 (cos(i;)) 3 sin(v) + 144 cos(v) sin(-u) + 768 sin(t>)(cos(i;)) 5 
-384 (cos(w)) 4 sin(v) + 288 (cos(v)) 2 sin(u) + 432 cos(u )« - 941 cos(v)t; 3 + 281 v 3 
+1344 (cos(v)) 5 t; + 1344 (cos(w)) 4 w - 1776 (cos(u )) 3 w - 768 (cos(v)) 6 v - 576 (cos(v)) 2 v 
+259 (cos(v)) V + 1481 (cos(v)) 3 v 3 

b\num = -24 - 800w(cos(w)) 3 sin(w) - 512 sin(w)w(cos(w)) 6 + 1088 w(cos(v)) 4 sin(u) 

-25 + 628(cos(v))V - 644 (cos(u))V + 216 u cos(w) sin(v) + 24wsin(w) 

+168 (cos(v)) 2 - 528 w(cos(w)) 2 sin(v) - 640 (cos(w)) V + 68 w 2 + 768 (cos(w)) 5 

+192 (cos('u)) 6 — 180 cos(f)w 2 + 385 (cos(v))V + 72 cos(«) - 456 (cos(v)) 3 - 

336 (cos(w)) 4 + 512 sin(w)w(cos(w)) 5 + 435 (cos(w)) 3 w 4 + 768 (cos(w))V 

-384(cos(w)) 7 + 192 v 2 (cos (v)) 7 - 192 v 2 (cos(v)) 6 - 75 cos(v)v 

b\ num = 24 + 1612f(cos(w)) 3 sin(w) + 736 sin(t> )w(cos(v)) 7 - 960 sin(w)w(cos(t;)) 6 

+676 w(cos(t;)) 4 sin(w) + 55 v 4 + 548 (cos(t> )) V + 804 (cos(u)) V 

— 240 v cos(u) sin(f ) — 32usin(t;) + 192 sin(u)(cos(v)) 8 w — 1050 (cos(v)) 4 v 4 

+672 (cos(w)) 8 - 300 (cos(w)) 2 + 136v(cos(t> )) 2 sin(v) - 1028 (cos(v)) V 

-60 v 2 + 132 (cos(w)) 5 - 1560 (cos(v)) 6 - 435 v 4 {cos(v)) 5 + 96 cos(w)w 2 

-256 (cos(v)) 8 v 2 - 265 (cos(w)) V - 24 cos(v) + 84 (cos(w)) 3 + 1164 (cos(w)) 4 

-2120 sin(vMcos(v)) 5 - 865 (cos(v)) 3 v 4 - 1264 (cos(v)) V - 64 (cos(v)) V 

-384 (cos(t>)) 7 + 192 (cos(w)) 9 + 448 v 2 (cos (v)) 7 + 776 w 2 (cos(t;)) 6 + 40 cos(v)v 4 
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b\ num = —72 — 3904u(cos(f )) 3 sin(v) — 1536 sva(v)v{cos(v)) 7 + 256 sin(u)i;(cos(u)) 6 
+512w(cos(w)) 4 sin(w) - 75 v 4 - 244(cos(w))V - 1980 (cos(v))V 
+648 w cos(w) sin(u) + 104 1; sin(u) + 1740 (cos(u)) 4 v 4 - 1536 (cos(v)) 8 + 792 (cos(v)) 2 
-944v(cos(w)) 2 sin(v) + 832 (cos(u))V + 156 v 2 + 960 (cos(t>)) 5 + 3648 (cos(w)) 6 
+580v 4 (cos(v)) 5 - 396 cos{v)v 2 + 512 (cos(v))V + 855 (cos(v))V + 120 cos(v) 
-696 (cos(v)) 3 - 2832 (cos(v)) 4 + 4864 sm{v)v{cos{v)f + 2165 (cos(u)) V 
+3168 (cos(iO)V - 384 (cos (v)) 7 - 192 v 2 (cos (v)) 7 - 1856w 2 (cos(w)) 6 - 225 cos(t> )v 
b\ nurn = 84 + 5416i;(cos(w)) 3 sin(v) + 640 sin(i;)t;(cos(i;)) 7 - 6720 sin(»i;(cos(») 6 
+4080 w(cos(v)) 4 sin(u) + 165 v 4 + 2540 (cos(w)) 3 w 2 + 2764 (cos(v)) 2 v 2 
-900wcos(w) sin(u) - 120vsin(t;) +2304 sin(?;)(cos(t;)) 8 v - 4720 (cos(v))V 
+1920(cos(w)) 8 - 1020 (cos(w)) 2 + 420 w(cos(w)) 2 sin(w) -920 w 4 (cos(f)) 6 
-4528 (cos(v)) V - 172 v 2 + 2736 (cos(v)) 5 - 4800 (cos(w)) 6 - 3380 w 4 (cos(t;)) 5 
+260 cos(v> 2 - 825 (cos(u)) V - 60 cos(t;) - 180 (cos(t>)) 3 + 3816 (cos(w)) 4 
-5120 sm(v)v(cos(v)) 5 - 3115 (cos(v)) 3 v 4 - 2912 (cos(v)) V - 768 (cos(v)) 9 v 2 
-4800 (cos(w)) 7 + 2304 (cos(v)) 9 + 2496 w 2 (cos (v)) 7 + 320 v 2 (cos(v)) 6 + 195 cos(t>)t> 4 

Since for small values of v, the above formulae are subject to heavy cancelations, we give the 
Taylor expansions of the coefficients b: 



_ 399187 52559 ^ 371082169 4 83360891 ^ 1467578899 8 

1 ~ 241920 ~ 304128 V + 58118860800 V ~ 523069747200 V ~ 355687428096000 V 
_ 17327 52559^ 2 3253170563 4 5070942803 ^ 86978398867 8 

2 ~ ~ 8640 + 38016 V ~ 14529715200 V + 261534873600 V ~ 88921857024000 V 
_ 597859 367913 2 2523373219 4 22329042629 6 1453392734357 8 

3 ~ 60480 ~ 76032 V + 2075673600 V ~ 130767436800 V + 88921857024000 V 

704183 367913 2 6122891963 4 133660742933 6 5024895032029 8 

4 ~ ~ 60480 + 38016 V ~ 2075673600 V + 261534873600 V ~ 88921857024000 V 

465133 1839565 , 3240803569 4 37612768013 fi 2927078073011 8 
6 5 = t> 2 H t; 4 H 

5 24192 152064 830269440 52306974720 35568742809600 

The PLTE of the method is given 

pUe = ( w 6 (6) + w 2 (10) ^2559_ {12) _52559_ 4 (8) \ 12 

^ V 912384 y 304128 * 912384 y 304128 y K ' 



3.5. Method PF-D3: Phase fitted + 1st, 2nd, 3rd Derivatives of Phase-lag are zero 

The fourth method of the family (PF-D3) is constructed by forcing the phase-lag function and 
its 1st, 2nd and 3rd derivatives to be zero at the frequency V = lo * h. The coefficients a are left 
the same, while for the coefficients b we have: 
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b = 



12 D 3 



h = 




where 



D 3 = v 5 ((cos(v)) 7 - (cos(V)) 6 - 3 (cos(v)f + 3 (cos(w)) 4 + 3 (cos(v)f - 3 (cos{v)) 2 
— cos(f) + 1) 



b i,num = -2096w 2 (cos(w)) 5 sin(t>) + 1508v 2 (cos(t> )) 4 sm(v) + 1648v 2 (cos(t;)) 3 sm(v) 

— 964 sin(-u)(cos(?;)) 2 + 36 f — 366 u 2 cos(f) sin(-u) + 135 v 5 cos(v) — 48 v 3 — 

384 (cos(v)) 8 f 3 — 336 (cos(i>)) 4 sin(i;) — 656 v 2 (cos(v)) 6 sin(-u) — 336 cos(v)v 3 

+1440 (cos(v)fv - 738 (cos{v )) 2 v - 384 sm(v){cos{v)) 7 + 168 (cos(v)) 2 sin(v) 

+234 cos(u)u + 528 (cos(u )) V - 1098 {cos{v)) 3 v - 2448 (cos(u))S 

+768 sm(v){cos{v)f + 192 sin(t>)(cos(t;)) 6 - 456 (cos(t;)) 3 sm(v) + 1008 {cos{v)) 3 v 3 

+72 cos(i;) sin(u) + 864t>(cos(i;)) 8 + 2286 (cos(/u)) 4 i> + 45 v 5 (cos(v)) 3 - 24 sm(v) 

+1200(cos(v))V - 1008 (cos(v)fv 3 + 45t> 5 + 135 v 5 (cos(v)) 2 

-576 (cos(t;)) 7 t; - 1296 (cos(u )) V + 832 sin(v)v 2 (cos(w)) 7 + 94 v 2 sin(v) 

+336 (cos(t;)) 7 t; 3 

ft |n«m = y - 692v 2 (cos(t;)) 5 sin(>) - 2407 u 2 (cos(t;)) 4 sin(u) + 235 t> 2 (cos(t;)) 3 sin(u) 

+776 -u 2 sin(i;)(cos(f)) 2 - 27 w + 123 v 2 cos(v) sin(u) - 90 cos(i>) + 48 1; 3 

-96 (cos(v))V + 852 (cos(v)) 4 sin(u) + 2560t; 2 (cos(t>)) 6 sin(t>) + 144 

-3690 (cos(u )) 5 f + 90 (cos{v)) 2 v - 96 sin(v)(cos(u )) 7 - 192 (cos(u)) 2 sin(v) — 144 cos{v)v 

-48 (cos(u)) V + 1458 (cos(u )) 3 v - 738 {cos{v)fv + 120 sin(w)(cos(u )) 5 

-1248 sin(w)(cos(t;)) 6 - 24 (cos(w)) 3 sm(v) - 912 sm(v)v 2 \cos{v)f - 720 (cos(v))V 

+360t;(cos(i;)) 8 + 315 {cos{v)) A v + 576 sin(w)(cos(t;)) 8 - 270 w 5 (cos(t>)) 3 - 1152 {cos{v)fv 

+12 sin(u) + 240 (cos(t;)) V + 1296 (cos(u)) V - 270v 5 (cos(v)) 2 - 90 (cos(v))V 

+3528 {cos{v)fv - 144 (cos(u)) 4 u 3 + 352 sm{v)v 2 {cos{v)) 7 + 288 (cos(v)) V 

-35 v 2 sin(v) - 1008 (cos(u))V 



and 
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b lnum = 3664v 2 (cos(V)) 5 sm(t;) + 3218 1; 2 (cos (t;)) 4 sin(u) - 812 v 2 (cos(v)) 3 sm(v) 
-1498 v 2 sin(i;)(cos(w)) 2 + 54 1> - 258 v 2 cos(v) sm(v) + 135 v 5 cos(d) - 84 v 3 
+ 1344 (cos 

(V))V _ 120 o (cos(v)) 4 sin(t;) - 2420 v 2 (cos(v)) 6 
+ 1344 sin(t>)?; 2 (cos(w)) 9 - 300 cos(v)v 3 + 4842 (cos(v)fv - 180 (cos(?;)) 2 i> 
+2496 sin(v)(cos(i;)) 7 - 384 (cos(v)) 10 v 3 + 312 (cos(w)) 2 sin(t;) + 324 cos(v)v 
+60 (cos(v)) 2 v 3 - 2502 (cos(v)) 3 v + 6498 (cos(v)fv - 1608 sin(t;)(cos(u)) 5 
+270 (cos(v)) 5 v 5 + 1488 sm(v)(cos(t;)) 6 + 264 (cos(v)) 3 sin(u) + 624 sm(t>)t> 2 (cos(t>)) 8 
+ 1092 (cos(v)) 3 v 3 - 6336v(cos(v)) 8 - 2052 (cos{v)) A v - 576 sin(t>)(cos(t;)) 8 
+855w 5 (cos(f)) 3 + 864 (cos(?;)) 9 t> - 24 sm(v) - 1644 (cos (v))V - 1476 (cos(v)fv 3 
+45 + 405v 5 (cos(t>)) 2 - 1152 sin(t;)(cos(t;)) 9 + 810 (cos(u))V - 3528 (cos(v)) 7 v 
+708 (cos(u)) 4 u 3 - 3920 sin(t;)t; 2 (cos(t;)) 7 - 192 (cos(v))V + 2016 v (cos (t>)) 10 
+58 v 2 sm(v) + 876 (cos(v)) 7 v 3 

b\ num = -5752w 2 (cos(t;)) 5 sin(t>) - 3789v 2 (cos(t> )) 4 sin(u) + 1333v 2 (cos(t>)) 3 sm(v) 

+2196 v 2 sin(v)(cos(w)) 2 - 81 u + 417v 2 cos(w) sin(t>) - 270 v 5 cos(v) + 108 v 3 

-2016(cos(u)) 8 u 3 + 1404 (cos(v)) 4 sin(t;) + 1036 t> 2 (cos(v)) 6 

-2048 sin(t>)v 2 (cos(f )) 9 + 468 cos(t>)t> 3 - 5310 (cos(v)fv + 234 (cos(v)) 2 v 

-4224 sin(t>)(cos(t;)) 7 + 576 (cos(v)) 10 v 3 - 432 (cos(u )) 2 sin(u) — 504 cos(v)v 

-36 (cos(u)) V + 3582 (cos(t>)) 3 t> - 10818 (cos(v)) 6 ^ + 2760 sin(i;)(cos(t;)) 5 

-1080 (cos(u ))V - 816 sin(v)(cos(t;)) 6 - 456 (cos(t;)) 3 sin(u) + 1360 sm(w)w 2 (cos(t;)) 

-1500 (cos(v)) 3 v 3 - 1152 (cos{v)) u v + 10224 t;(cos(t;)) 8 + 3609 (cos(v)) 4 v 

-960 sin(t;)(cos(w)) 8 - 1170 w 5 (cos(t;)) 3 + 2592 (cos(t;)) 9 w + 36 sm(v) 

+768 sin(t;)(cos(i;)) 10 + 2484 (cos(t;)) 6 t; 3 + 1500 (cos(v)) 5 v 3 + 192 (cos(v)) ll v 3 

-810w 5 (cos(w)) 2 - 704 sin(v)« 2 (cos(v)) 10 + 1920 sin(t;)(cos(i;)) 9 - 1350 (cos(v))V 

+792 (cos(t;)) 7 t; - 1116 (cos(u)) 4 u 3 + 6032 sin(t;)t; 2 (cos(t;)) 7 - 360 v 5 (cos(v)f 

-480 (cos(u)) V - 3168 v (cos(w)) 10 - 81 v 2 sin(v) - 180 (cos(u)) V 
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bl num = 13216 v 2 {cos{v)fs\n{v) + 10404 v 2 (cos (v) ) 4 sm(v) - 4672 v 2 (cos(v )) 3 sin(v) 

-4932 sin(t;)(cos(t;)) 2 + 180 u - 762 v 2 cos(w) sin(t>) + 405v 5 cos(v) - 240 v 3 

+3456 (cos(u))V - 3888 (cos(v )) 4 sin(u) - 4048v 2 (cos(w)) 6 sin(v) 

+2816 sin(v)t; 2 (cos(t;)) 9 - 912 cos{v)v 3 + 15192 (cos(w)) 5 t> + 162 (cos(v)) 2 v 

+7680 sin(t;)(cos(t>)) 7 - 768 (cos(v)) 10 v 3 + 1080 (cos(u)) 2 sin(v) + 1062 cos(t;)t; 

-432 (cos(u )) V - 9054 (cos(t>)) 3 t> + 22680 (cos(t;)) 6 t; - 5856 sin(?;)(cos(t;)) 5 

+2160 (cos(v)) V + 2112 sin(t;)(cos(t;)) 6 + 1320 (cos(t;)) 3 sm(v) - 4288 sin(v)t> 2 (cos(t;)) 8 

+3504 (cos(v)) 3 v 3 + 4608 (cos(v)) u i; - 72 cos(f) sin(v) — 17856 v(cos(v)) 

-9774 (cos(t;)) 4 t; + 3840 sin(t>)(cos(t;)) 8 + 3375 (cos (t> )) 3 - 9216 (cos(v)) 9 v 

-72 sin(v) - 3072 sin(t;)(cos(t>)) 10 - 5520 (cos(v)fv 3 - 4272 (cos(v))V + 135 v 5 

-768 (cos(v)) n u 3 + 1485t> 5 (cos(V)) 2 + 2816 sin(t>)?; 2 (cos(f )) 10 

-3072 sin(t;)(cos(7;)) 9 + 3600 (cos(v))V - 2592 (cos(v)) 7 v + 3504 (cos(u)) 4 u 3 

-10688 sin(t;)t; 2 (cos(t;)) 7 + 360v 5 (cos(t;)) 7 + 1080v 5 (cos(w)) 6 

+1536 (cos(u)) 9 u 3 + 4608 v(cos(v)) w + 138 v 2 sin(v) + 912 (cos(u))V 

Since for small values of v, the above formulae are subject to heavy cancelations, we give the 
Taylor expansions of the coefficients b: 



_ 399187 52559 2 11315653 4 5807033 6 614853845 

1 ~ 241920 ~ 228096 V + 1937295360 V ~ 13076743680 V ~ 17072996548608 
17327 52559 2 5758537 4 225159101 6 8200289261 8 

2 " ~ 8640 + 28512 V ~ 14676480 V + 6538371840 V ~ 4268249137152 V 
_ 597859 _ 367913 ^ 154801723 ^ _ 2799488011 ^ 1063054198007 ^ 8 

3 ~ 60480 ~ 57024 V + 69189120 V ~ 6538371840 V + 21341245685760 V 
704183 367913 2 381346481 4 9217976399 6 5100295346143 

4 ~ ~ 60480 + 28512 V 69189120 V + 6538371840 V ~ 21341245685760 1 
_ 465133 1839565 2 67543471 4 241481599 6 16316044646989 8 

5 ~ 24192 114048 V + 9225216 V ~ 118879488 V + 42682491371520 V 



The PLTE of the method is given 

( 



/ 52559 o /-irA 52559 « /a\ 52559 a (q\ 

plte = w 2 y {10 > H w s y {4 > H wV 8) 

1 228096 y 912384 y 152064 y 



52559 ^6 (6) , 52559 (12) \ h i2 
228096 y 912384 y 



3.6. Method PF-D4: Phase fitted + 1st, 2nd, 3rd, 4th Derivatives of Phase-lag are 

zero 

The fifth method of the family (PF-D4) is constructed by forcing the phase-lag function and 
its 1st, 2nd, 3rd and 4th derivatives to be zero at the frequency V = u> * h. The coefficients a are 
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left the same, while for the coefficients b we have: 



bo 
h 





b 2 = 




where 



D 4 = v 6 (cos (v) + 1) (sin (v)) 5 



and 



b\ num = 3069 sin(u)(cos(v))V - 6408 sin(v)(cos(v))V - 4218 sin(v)(cos(v)) V 

—864 sin(v) cos(-u)i> 4 + 1728 sin(u)(cos(f )) 3 v A + 2016 sin(v)(cos(v)) 4 i; 4 

-864 sin(f)(cos(w)) 5 w 4 + 1380 sin(v ) cos(v)v 2 - 1152 sin(u)(cos(t> )) 2 v 4 

+2568 sin(v)(cos(w)) 5 w 2 + 3408 sin(v)(cos(w))V - 960 sin(v)(cos(w))V 

-144 v + 1440 (cos(w)) 4 sin(w) + 1029 cos(v)t; 3 + 60 sin(u) + 372 v 3 + 4836 (cos(v)) 4 v 3 

-5952 (cos(v)fv + 1728 (cos(v)) 2 v - 540 (cos(w)) 2 sin(v) — 648 cos(v)v — 3234 (cos(u)) v 

+3912 (cos(v)) 3 v + 1728 (cos(t;)) 6 t; - 480 sin(w)(cos(w)) 5 + 2688 (cos(v)) 7 v 

-2064 (cos(v)) 6 u 3 + 600 (cos(v)) 3 sin(u) - 204 sm(v)v 2 - 4786 (cos(v)) 3 v 3 

+96 sin(w)t; 4 - 120 cos(w) sm(v) - 3312 (cos(w)) 4 t; + 6176 (cos(v)) 5 v 3 

-960 sin(w)(cos(w)) 6 - 2464 (cos(«))V 

b\ nurn = -1047 sin(u)(cos(v))V + 3552 sin(u)(cos(v)) V - 2295 sin(w)(cos(w)) 3 w 2 

+48 sin(v) cos(f )f 4 + 384 sin(v)(cos(w)) 3 u 4 — 912 sin(f )(cos(w)) 4 v 4 

—912 sin(f )(cos(v)) 5 u 4 + 78 sin(v) cos(v )v 2 + 384 sin(f ){cos{v)) 2 v A 

+5184 sin(v)(cos(v))V - 2208 s\n{v)(cos{v)fv 2 - 2832 sm{v)v 2 (cos{v)) 7 

+480 sin(v)(cos(v))V + 12 u - 600 (cos(w)) 4 siii(v) - 131 cos{v)v 3 - 35 v 3 

+2929 (cos(w)) 4 w 3 + 3096 (cos(w)) 5 w + 528 (cos{v)) 2 v + 120 (cos(w)) 2 sm(v) + 120 cos(v)u 

-433 (cos(v)) 2 v 3 - 1584 (cos(w)) 3 w + 5472 (cos(v)fv - 1440 sin(w)(cos(w)) 5 

-1632 (cos(w)) 7 w - 4128 (cos(v)fv 3 + 540 (cos(w)) 3 sin(v) - 27 sin(v)v 2 

+2181 (cos(v))V + 48 sm(v)v A - 60 cos{v) sm(v) - 3516 (cos(w)) 4 w - 3480 (cos(v)) 5 u 3 

+480 sin(w)(cos(w)) 6 + 1520 (cos(v)) V + 1712 (cos(v)) 8 v 3 

-2496 (cos(v)) s v + 960 sin(t>)(cos(w)) 7 + 480 siii(v)(cos(«)) V 

6| num = 2997 sin(u)(cos(v))V + 3906 sin(v)(cos(u))V + 498 sin(v)(cos(u))V 

—936 sin(v) cos(v)v 4 — 48 sin(f )(cos(v)) 3 v 4 + 24 sin(u)(cos(f )) 4 t> 4 

+2904 sin(f)(cos(w)) 5 w 4 + 1632 sin(w) cos(w)w 2 - 1008 sin(v)(cos(u)) V 

-15396 sin(v)(cos(u))V - 21792 sin(v)(cos(u)) V - 7488 (cos(v))V 
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+11376 sin(V)t> 2 (cos(t>)) 7 + 2880 sin(v)(cos(v)) V - lUv - 1800 (cos(V)) 4 sin(u) 

+1101 cos(t;)t; 3 + 60 sin(v) + 300 v 3 - 4740 {cos(v)) 4 v 3 + 12768 (cos{v)) 5 v 

+1152 (cos(v)) 2 v - 180 (cos(v)) 2 sin(w) — 792 cos(t> )v — 2910 (cos(i> )) v 

+1464 (cos(t;)) 3 t; - 15552 (cos(v)fv + 3120 sin(t;)(cos(y)) 5 - 27264 (cos(v)) 7 v 

+13632 (cos(v)) 6 t; 3 - 120 (cos(i;)) 3 sin(t;) - 168 sm(v)v 2 - 3244 (cos(v)) 3 v 3 

+24 sin(v)v 4 - 120 cos(u) sin(u) + 5328 (cos(u )) 4 v - 5980 (cos(v)fv 3 

+7680 sin(t;)(cos(V)) 6 + 15296 (cos(v))V - 6912 (cos(V)) 8 t> 3 + 9216 (cos(v)) 8 t; 

-2880 sin(v)(cos(?;)) 7 + 14112 sin(u)(cos(t;)) 8 t; 2 + 13824 (cos {v)) 9 v 

-1920 sin(v)(cos(u))V - 1920 sin(u)(cos(u ))V - 5760 sin(t>)(cos(t;)) 8 

b\ num = -3321 sin(u)(cos(v))V + 8268 sin(u)(cos(?;))V - 6861 sin(u)(cos(u))V 

+108 sin(v) cos(?;)t> 4 + 1224 sin(-u)(cos(u)) V - 1812 sm(tO(cos(t;)) V 

-1812 sin(u)(cos(V)) V + 162 sin(u) cos(v)v 2 + 1224 sin(v)(cos(/y)) V 

+10188 siii(v)(cos(v))V + 3552 sin(v)(cos(w))V + 960 sin(v)(cos(w)) 9 t; 4 

+3648 (cos(w))V + 5328 s\n(v)v 2 {cos{v)) 7 - 7872 sin(w)w 2 (cos(i;)) 9 

-480 sm.(v)(cos{v)fv 4 + mv - 1320 (cos(w)) 4 sin(u) - 297 cos(w)w 3 - 81 u 3 

+7815 (cos(w)) 4 w 3 + 4104 (cos(w))S + 1728 (cos(w)) 2 w + 360 (cos(u)) 2 sin(u) 

+360 cos(u )v - 1287 (cos(v)) V - 4464 {cos{v)) 3 v + 5520 (cos(u )fv 

-2160 sin(t;)(cos(i;)) 5 + 5760 (cos(v)) 7 v - 5828 (cos(u))V + 1380 (cos(i;)) 3 sin(u) 

-81 sm(v)v 2 + 6315 {cos{v)) 3 v 3 + 108 sm{v)v 4 - 180 cos(t;) sin(v) - 9396 {cos{v)) 4 v 

-5244 (cos(v))V - 960 sin(t>)(cos(t;)) 6 - 3792 (cos(u))V + 3840 sin(w)(cos(v)) 9 

+3904 (cos(v)) 10 v 3 - 4208 (cos(t;)) 8 t; 3 + 10560 (cos{v)fv - 2880 sin(t;)(cos(w)) 7 

-6528 $m(v)(cos(v)fv 2 - 5760 (cos(t;)) 9 t; - 480 sm(v)(cos(v)) 7 v 4 

+960 sin(i;)(cos(t;)) V + 1920 sin(t;)(cos(i;)) 8 - 8448 (cos(t;)) 10 u 
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b\ num = 9207 sm{v)(cos{v)) 2 v 2 + 21552 sin(v)(cos(v))V + 6498 sm{v)(cos{v)fv 2 

-2784 sin(t>) cos(v)t; 4 - 2112 sin(u)(cos(v)) 3 u 4 - 1824 sin(V)(cos(v)) V 

+11040 sin(v)(cos(u))V + 5148 sin(u) cos(v)v 2 - 3072 sin(v)(cos(u))V 

-60312 sin(v)(cos(u))V - 81288 sin(u)(cos(v)) V - 1536 sin(v)(cos(u))V 

-14592 (cos(v))V + 27888 sin(?;)t> 2 (cos(f )) 7 + 11328 sin(v)v 2 (cos(t;)) 9 

+10944 sin(u)(cos(w))V - 432 u - 8160 (cos(w)) 4 sin(v) + 3183 cos(v)v 3 

+180 sin(v) + 828 v 3 - 23348 (cos(v))V + 55680 (cos(v)) 5 v + 2880 {cos(v)) 2 v 

-180 (cos(w)) 2 sin(v) - 2520 cos(v)v - 8598 (cos(v)) V + 1560 (cos(v)) 3 v 

-54720 (cos(w)) 6 w + 12000 sin(w)(cos(w)) 5 - 93120 (cos(v)) 7 v + 47712 (cos(v))V 

-1080 (cos(u)) 3 sin(v) - 468 sin(t> )v 2 — 7742 (cos(f )) 3 u 3 + 96 sin(t;)i; 4 

-360 cos(v) sin(v) + 24240 (cos(w)) 4 w - 29032 (cos(u)) V + 27360 sin(w)(cos(w)) 6 

+53008 (cos(v)) 7 v 3 - 3840 sin(w)(cos(w)) 9 - 6016 (cos(v)) w v 3 

-13728 (cos(v)) 8 v 3 + 17280 (cos(v)) s v - 6720 sin(w)(cos(w)) 7 

+32832 sin(v)(cos(u))V + 23040 (cos(v)) 9 v - 6400 (cos(u)) n t; 3 

+13440 sin(v)(cos(u)) 1( V - 4608 sin(u)(cos(w)) V - 4608 sin(v)(cos(u)) 8 u 4 

-11520 sin(w)(cos(w)) 8 - 7680 sin(w)(cos(w)) 10 + 10752 (cos(v)) 10 v 

-1536 sin(v)(cos(v)) 10 v 4 + 15360 (cos(v)) n u 

Since for small values of v, the above formulae are subject to heavy cancelations, we give the 
Taylor expansions of the coefficients b: 



399187 262795 , 17265277 4 38566679 e 52935007231 8 

bi = v 2 H w 4 v 6 v 8 

241920 912384 4358914560 38041436160 426824913715200 
_ 17327 262795 ^ 2 2649128441 ^ 4 2650726483 ^ 374485131133 ^ 

2 ~ ~ 8640 + 114048 V ~ 4358914560 V + 52306974720 V ~ 106706228428800 V 
_ 597859 1839565 2 1110676079 4 6919361527 6 1637603830619 8 

3 ~ 60480 ~ 228096 V + 311351040 V ~ 8047226880 V + 15243746918400 V 
704183 1839565 2 5518849841 4 14263211663 6 73005517242211 

4 ~ ~ 60480 + 114048 V 622702080 V + 4755179520 V ~ 106706228428800 
_ 465133 9197825 2 734695627 4 183227481067 6 9908801489731 8 

5 ~ 24192 ~ 456192 V + 62270208 V 41845579776 V + 8536498274304 V 



if 



The PLTE of the method is given 



. 52559 in m 262795 fi m 

plte = ( w 10 v {2 > H w 6 y {6 > + 

912384 456192 
262795 o m 262795 4 f8 > 52559 262795 2 rin v ll2 

w 8 y {4 > + u>V 8) + y (12) + w 2 y (loy )h 12 

912384 y 456192 y 912384 y 912384 tf ; 
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4. Stability Analysis 

The stability of the new methods is studied by considering the test equation 

d 2 y(t) 



dt 2 



-a 2 y(t) (23) 



and the linear multistep method ([3]) for the numerical solution. In the above equation a ^ u (u 
is the frequency at which the phase-lag function and its derivatives vanish). Setting s = ah and 
v = u>h we get for the characteristic equation of the applied method 

J/2 

£ Ms^vW + z-i) + A (s\v) = (24) 
i=i 

where 

Aj{s 2 ,v) = aj_ J (v) + s 2 - b j_.(v) (25) 

The motivation of the above analysis is straightforward: Although the coefficients of the method 
(|3|) are designed in a way that the phase-lag and its first derivatives vanish in the frequency u, the 
frequency u itself is unknown and only an estimation can be made. Thus, if the correct frequency 
of the problem is a we have to check if the method is stable, that is if the roots of the characteristic 
equation lie on the unit disk. For this reason we draw at the v-s plane the areas in which the 
method is stable. Figure [T] shows the stability region for the six methods (the classical one, the 
phase fitted one and those with first, second, third and fourth phase lag derivative elimination). 
Note here that the s-axis corresponds to the real frequency while the u-axis corresponds to the 
estimated frequency used to construct the parameters of the method. 



5. Numerical Results 

Numerical experiments have been carried out for two orbital problems. Since the classical 
method is well studied, we only present the new methods in comparison to the classical one. 



5.1. The 2-Body Problem 

In this problem we test the motion of two bodies in a reference system that is fixed in one of 
them. Moreover, the motion is planar, thus, we only have to calculate the x and y coordinates of 
the second body. The differential equations are: 

X = -(3,2+^2)3 

•■ y 

y — JxZ+yZf 
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and the initial conditions are 

x(0) = 1 - e ir(0) = 



j/(0) =0 y(0) = 




where e is the eccentricity. Figure [2] presents the accuracy of the methods expressed by —logio 
(error at the end point which is after 100 periods) versus the logio (total steps). The results clearly 
show the improvement in the accuracy using the proposed methodology. 



5.2. The 5-Outer Planet System 

The next problem concerns the motion of the five outer planets relative to the sun. The problem 
falls in the category of the N-Body problem which is the problem that regards the movement of N 
bodies under Newton's law of gravity. It is expressed by a system of vector differential equations 

y> = G £ (26) 

where G is the gravitational constant, rrij is the mass of body j and yi is the vector of the position 
of body i. It is easy to see that each vector differential equation of ([26]) can be analyzed into three 
simplified differential equations, that express the three directions x, y, z. 

The above system of ODEs cannot be solved analytically. Instead we produce a highly accurate 
numerical solution by using a 10-stage implicit Runge-Kutta method of Gauss with 20th algebraic 
order, that is also symplectic and A-stable. The method can be easily reproduced using simplifying 



assumptions for the order conditions (see iButcherl (|2003l )). The reference solution is obtained 



by using the previous me thod to integrate t he N-body problem for a specific time-span and for 



different step-lengths. In lHairer et al.l ([2002) the data for the five outer planet problem is given 
(these data are summarized in table [1]). Masses are relative to the sun, so that the sun has mass 
1. In the computations the sun with the four inner planets are considered one body, so the mass 
is larger than one. Distances are in astronomical units, time is in earth days and the gravitational 
constant is G = 2.95912208286 • 10" 4 . The system of equations (J26]) has been solved for t G [0, 10 6 ], 
for which time-span, the previously mentioned method of Gauss produces a 10.5 decimal digits 
solution. Figure [3] presents the accuracy of the methods expressed by —logio (error at the end 
point) versus the logio (total steps). The results clearly show the improvement in the accuracy 
using the proposed methodology. 



5.3. The figures 



In Figure [T] we see the stability regions (v-s plane) for the classical Quinlan-Tremaine method, 
the PF-D0, PF-D1, PF-D2, PF-D3 and PF-D4 methods (from left to right and from top to bottom). 
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In Figure [2] we see the accuracy of the methods expressed by —logio (error at the end point which 
is after 100 periods) versus the logio (total steps) in the two body problem for eccentricity e = 0.5 
and step size h = 0.1. In Figure [3] we see the accuracy of the methods expressed by —logio (error at 
the end point which is after 10 6 days) versus the logio (total steps) for the 5-Outer Planet System. 



s-v plane QT10-PF-DO (v— S plane) QT10-PF-D1 (v-s plane) 




0.0 0.5 1.0 IS 2.0 2.5 3.0 0.0 0.5 1.0 1.5 2.0 2,5 3.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 



Fig. 1. — Figure 1. 



6. Conclusions 



We have presented a new family of 10-step symmetric multistep numerical methods with 
improved characteristics concerning orbital problems. The methods were constructed by adopting a 
new methodology which, apart from the phase fitting at a predefined frequency, it eliminates the first 
derivatives of the phase lag function at the same frequency. The result is that the phase lag function 
becomes less sensitive on the frequency near the predefined one. This behavior compensates the fact 
that the exact frequency can only be estimated. Experimental results demonstrate this behavior 
by showing that the accuracy is increased as the number of the derivatives that are eliminated is 
increased. 
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Fig. 2. — Figure 2. 
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Table 1: Initial Data for the 5-outer planet problem 
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